function [] = computeElasticities_consumerAdoptions(runId)

%%% Add paths
addpath('../model_funcs');
addpath('../general_funcs');

%%% Get input / output paths
config = getConfig_consumerAdoptions(runId);

if config.modelId ~= 3
	error('This code is only for modelId == 3 (continuous time models).');
end

outputFolder = createFolderIfNotExist(sprintf('%s/elasticities', config.runFolder));
outputMatFile = sprintf('%s/elasticities.mat', outputFolder);
outputCSVFile1  = sprintf('%s/elasticities.csv', outputFolder);
outputCSVFile2  = sprintf('%s/elasticitiesStrings.csv', outputFolder);

%%% Load data
load(sprintf('%s/runData.mat', config.runFolder), 'data');
%% Get paramNames
paramNames = get_param_names_consumerAdoptions(data);

%%%  Load parameter estimates from the 3 submodels
load(sprintf('%s/results.mat', config.runFolder), 'params_hat', 'params_ses');
params = params_hat;
%%% Split params_parts + save everything in object paramsParts
paramsParts  = split_combine_parts(data.dims,  1, params);


%%% Read dimensions
NumSpells = data.dims.dims1{4};
K1 = data.dims.dims1{5};
K2 = data.dims.dims1{6};
NumParams = data.dims.NumParams;
NumParts = length(data.Xparts);
	
myfunc1 = @(x) x; % Identity function
myfunc2 = @(x) 1; % one
myfunc3 = @(x) 1-exp(-x); % log(1+x) function

elastNames = {'wrt_X', 'wrt_eX', 'wrt_eX_minus1', 'wrt_Area_TIMES_eX_minus1'};


%%%% IDEA: for each spell, compute average elasticities across event types (adoptions in canton k), then take a weighted average across
% spells, where the weight is the duration of the spell.

elasticities1 = zeros(NumParams, NumSpells);
elasticities2 = zeros(NumParams, NumSpells);
elasticities3 = zeros(NumParams, NumSpells);

% Assign things that won't depend on spell
data_uu.surfaces                 = data.surfaces;
data_uu.dims.dimType             = data.dims.dimType;
data_uu.dims.NumParts            = data.dims.NumParts;
data_uu.dims.NumParams           = data.dims.NumParams;
data_uu.dims.Xpart_2_NumX        = data.dims.Xpart_2_NumX;
data_uu.dims.Xpart_2_NumX_FEs    = data.dims.Xpart_2_NumX_FEs;
data_uu.dims.Xpart_2_Num_FE_vals = data.dims.Xpart_2_Num_FE_vals;
data_uu.dims.NumFEvals2Keep      = data.dims.NumFEvals2Keep;
data_uu.dims.NumObs   = K1*K2;
data_uu.dims.dims1{1} = K1*K2;
data_uu.dims.dims1{2} = K1;
data_uu.dims.dims1{3} = K2;
data_uu.dims.dims1{4} = 1;
data_uu.dims.dims1{5} = K1;
data_uu.dims.dims1{6} = K2;

for pp = 1:NumParts
	data_uu.Xparts{pp}.NumX_FE_vals = data.Xparts{pp}.NumX_FE_vals;
	data_uu.Xparts{pp}.Xnames       = data.Xparts{pp}.Xnames;
	data_uu.Xparts{pp}.X_FE_labels  = data.Xparts{pp}.X_FE_labels;
end
data_uu.Xparts{1} = data.Xparts{1};
data_uu.Xparts{5} = data.Xparts{5};
data_uu.Xparts{6} = data.Xparts{6};

X2 = data.Xparts{2}.X; % (NumSpells*K1) x NumX2
X2 = reshape(X2, [NumSpells K1 size(X2,2)]); % NumSpells x K1 x NumX2
X2 = permute(X2, [2 3 1]); % K1 x NumX2 x NumSpells

X3 = data.Xparts{3}.X; % (NumSpells*K2) x NumX3
X3 = reshape(X3, [NumSpells K2 size(X3,2)]); % NumSpells x K2 x NumX3
X3 = permute(X3, [2 3 1]); % K2 x NumX3 x NumSpells

X_FEs2 = data.Xparts{2}.X_FEs; % (NumSpells*K1) x NumX_FEs2
X_FEs2 = reshape(X_FEs2, [NumSpells K1 size(X_FEs2,2)]); % NumSpells x K1 x NumX_FEs2
X_FEs2 = permute(X_FEs2, [2 3 1]); % K1 x NumX_FEs2 x NumSpells

X_FEs3 = data.Xparts{3}.X_FEs; % (NumSpells*K2) x NumX_FEs3
X_FEs3 = reshape(X_FEs3, [NumSpells K2 size(X_FEs3,2)]); % NumSpells x K2 x NumX_FEs3
X_FEs3 = permute(X_FEs3, [2 3 1]); % K2 x NumX_FEs3 x NumSpells

denom = data.surfaces; % K x 1

tic
for uu = 1:NumSpells
	data_uu.Xparts{2}.X     = X2(:,:,uu);     % K1 x NumX
	data_uu.Xparts{2}.X_FEs = X_FEs2(:,:,uu); % K1 x NumX_FEs
	
	data_uu.Xparts{3}.X     = X3(:,:,uu);     % K2 x NumX
	data_uu.Xparts{3}.X_FEs = X_FEs3(:,:,uu); % K2 x NumX_FEs
	
	data_uu.Xparts{4}.X     = data.Xparts{4}.X(uu,:); % 1 x NumX
	data_uu.Xparts{4}.X_FEs = data.Xparts{4}.X_FEs(uu,:); % 1 x NumX_FEs

	elasts1_uu = compute_mean_Xfunc_times_beta(data_uu.dims, data_uu.Xparts, myfunc1, paramsParts)';
	elasts2_uu = compute_mean_Xfunc_times_beta(data_uu.dims, data_uu.Xparts, myfunc2, paramsParts)';
	elasts3_uu = compute_mean_Xfunc_times_beta(data_uu.dims, data_uu.Xparts, myfunc3, paramsParts)';

	% Compute elasts4 wrt Area_TIMES_eX_minus1 (for variables X = log([1+Xtilda]/Area) )
	Xfunc_mean_beta_parts = cell(NumParts,1);
	for pp = 1:NumParts
		Xpart_pp = data_uu.Xparts{pp};
		if pp == 1 || pp == 4 % K1-K2 part, T-part
			Xfunc_mean_beta_parts{pp}.X_Y = NaN * ones(data_uu.dims.Xpart_2_NumX(pp),1); % NumX x 1
		end
		if pp == 2 % T-K1 part
			NumXp = size(Xpart_pp.X, 2);
			X_pp = reshape(Xpart_pp.X, [1 K1 NumXp]);
			Xfunc_mean_beta_parts{pp}.X_Y = paramsParts{pp}.beta .* mean(reshape(1 - exp(-X_pp)./denom', [K1 NumXp]), 1)'; % NumXp x 1
		end
		if pp == 3 % T-K2 part
			NumXp = size(Xpart_pp.X, 2);
			X_pp = reshape(Xpart_pp.X, [1 K2 NumXp]);
			Xfunc_mean_beta_parts{pp}.X_Y = paramsParts{pp}.beta .* mean(reshape(1 - exp(-X_pp)./denom', [K2 NumXp]), 1)'; % NumXp x 1
		end
		if pp == 5 % K1-part
			Xfunc_mean_beta_parts{pp}.X_Y = paramsParts{pp}.beta .* mean(1 - exp(-Xpart_pp.X)./denom, 1)'; % NumX x 1 
		end
		if pp == 6 % K2-part
			Xfunc_mean_beta_parts{pp}.X_Y = paramsParts{pp}.beta .* mean(1 - exp(-Xpart_pp.X)./denom, 1)'; % NumX x 1 
		end
		NumX_FEs = size(Xpart_pp.X_FEs,2);
		Xfunc_mean_beta_parts{pp}.X_FEs_Y = cell(NumX_FEs,1);
		for ff = 1:NumX_FEs
			NumFE_vals_ff = Xpart_pp.NumX_FE_vals(ff);
			Xfunc_mean_beta_parts{pp}.X_FEs_Y{ff} = NaN * ones(NumFE_vals_ff,1);
		end
	end
	elasts4_uu = split_combine_parts(data_uu.dims, 3, Xfunc_mean_beta_parts)'; % 1 x NumParams
	
	% Store everything
	elasticities1(:,uu) = elasts1_uu;
	elasticities2(:,uu) = elasts2_uu;
	elasticities3(:,uu) = elasts3_uu;
	elasticities4(:,uu) = elasts4_uu;

	displayRemainingTime(uu, NumSpells, 100);
end

spellWeights = data.spellDurations/sum(data.spellDurations); % NumSpells x 1
elasts1 = elasticities1*spellWeights;
elasts2 = elasticities2*spellWeights;
elasts3 = elasticities3*spellWeights;
elasts4 = elasticities4*spellWeights;

allElasticities = [elasts1 elasts2 elasts3 elasts4]; % NumParams x 4

% Compute standard errors on elasticities
allElasticities_ses = allElasticities .* params_ses./params;


%%%% Get the relevant elasticity for each variable, depending on any transformation applied
paramNames = strrep(paramNames, 'Beta ', '');
paramNames = strrep(paramNames, ' on NewConsumers', '');
rowNamesMapping = getRowNamesMapping();
[flag,idxes] = ismember(paramNames, rowNamesMapping(:,1));
idxes                = idxes(flag);
allElasticities      = allElasticities(flag,:);
allElasticities_ses  = allElasticities_ses(flag,:);
paramNames           = rowNamesMapping(idxes,2);

myVariablesAndElastType = {
	'Log density of providers at destination at k2'          , 2;
	'Log density of consumers living in origin at k1'        , 4;
	'Log density of public transit stops at k2'              , 4;
	'Log density of rental agencies at k2'                   , 4;
	'Log surface area of destination canton at k2'           , 2;
	'Income (in 100,000 euros) at k1'                        , 1;
	'Unemployment rate at k1'                                , 1;
	'Frac. of people without high school diploma at k1'      , 1;
	'Frac. of people with higher educ. diploma at k1'        , 1;
	'Fraction of votes for ecological candidate at k1'       , 1; 
	'Google Trend of competitor at k1'                       , 1; 
	'Precipitation at destination (in decimeters) at k2'     , 1;
	'Temperature at destination (in celsius) at k2'          , 1;
};
[flag2,idxes] = ismember(myVariablesAndElastType(:,1), paramNames); assert(all(flag2));
paramNames = paramNames(idxes);
allElasticities     = allElasticities(idxes,:);
allElasticities_ses = allElasticities_ses(idxes,:);
NumElasts = length(myVariablesAndElastType);
myelasts = zeros(NumElasts,1);
myelasts_ses = zeros(NumElasts,1);
for ee = 1:NumElasts
	myelasts(ee)     = allElasticities(ee,myVariablesAndElastType{ee,2});
	myelasts_ses(ee) = allElasticities_ses(ee,myVariablesAndElastType{ee,2});
end

% Rename those paramNames with log-density (after the fact)
paramNames = strrep(paramNames, 'Log density', 'Density');
paramNames = strrep(paramNames, 'Log surface', 'Surface');

M1 = table(myelasts, myelasts_ses);
M1.Properties.RowNames = paramNames;
writetable(M1, outputCSVFile1, 'WriteRowNames', true);

M2 = cell(NumElasts, 1);
for ee = 1:NumElasts
	M2{ee} = sprintf('%.2f (%.2f)', myelasts(ee), myelasts_ses(ee));
end
M2 = cell2table(M2);
M2.Properties.VariableNames = {'Elasticity'};
M2.Properties.RowNames = paramNames;
disp(M2);
writetable(M2, outputCSVFile2, 'WriteRowNames', true);


%%%%% Save to output file
save(outputMatFile, 'M1', 'M2');

end
